rsf_2023_working.R

Resource selection functions

Björn Reineking, 2023-08-05

Model the relative density of animals (also called range distribution or utilisation distribution) as a function of environmental predictors.

We will use the buffalo data set. Install animove metapackage https://github.com/AniMoveCourse/animove_R_package install.packages(“remotes”) remotes::install_github(“AniMoveCourse/animove_R_package”)

This code needs ctmm version 1.1.1 (or newer) remotes::install_github(“ctmm-initiative/ctmm”) Loading packages

library(animove)
library(ctmm)
library(sf)
library(mgcv)
library(mvtnorm)
library(lubridate)

Load buffalo data

See Kami’s slides for how we got here…

data(buffalo_utm)

Environmental data: topography, waterways, and NDVI

data(buffalo_env)
raster::plot(buffalo_env)

Animals and elevation

raster::plot(raster(buffalo_env, 1))
points(buffalo_utm)

## Loading required package: move

## Loading required package: geosphere

## Loading required package: sp

## Loading required package: raster

## 
## Attaching package: 'raster'

## The following object is masked from 'package:nlme':
## 
##     getData

## The following objects are masked from 'package:ctmm':
## 
##     projection, projection<-

## 
## Attaching package: 'move'

## The following object is masked from 'package:ctmm':
## 
##     speed

To speed up analyses, we will only work with one individual, Cilla

cilla <- buffalo_utm[["Cilla"]]
raster::plot(raster(buffalo_env, 1), ext = extent(cilla) * 2)
lines(cilla)

The first day of Cilla shows some unrealistic movements. We remove those observations

cilla <- cilla[timestamps(cilla) > min(timestamps(cilla)) + days(1), ]

raster::plot(raster(buffalo_env, 1), ext = extent(cilla) * 2)
lines(cilla)

Minimal example of rsf.fit

Create telemetry object

cilla_telemetry <- as.telemetry(cilla)

## Minimum sampling interval of 3 minutes in Cilla

Fit ctmm model

cilla_guess <- ctmm.guess(cilla_telemetry, CTMM=ctmm(isotropic = TRUE), interactive = FALSE)
cilla_select <- ctmm.select(cilla_telemetry, cilla_guess)

Fit akde

cilla_akde <- akde(cilla_telemetry, cilla_select)

## Default grid size of 3 minutes chosen for bandwidth(...,fast=TRUE).

plot(cilla_akde)

Create named list of rasters

be <- list("elev" = raster(buffalo_env, "elev"),
           "slope" = raster(buffalo_env, "slope"),
           "var_NDVI" = raster(buffalo_env, "var_NDVI"))

The integrator = “Riemann” option is still in testing, we use it here because it is much faster

cilla_rsf_riemann <- rsf.fit(cilla_telemetry, cilla_akde, R = be, integrator = "Riemann")

## Warning in .local(x, ...): This function is only useful for Raster* objects
## with a longitude/latitude coordinates

## Maximizing likelihood @ n=166035

## Calculating Hessian

summary(cilla_rsf_riemann)

## $name
## [1] "OUF"
## 
## $DOF
##      mean      area diffusion     speed 
##  8.323012 16.577960 72.507827 71.911902 
## 
## $CI
##                                             low          est         high
## var_NDVI (1/var_NDVI)             -188.04790701 -72.68967041  42.66856618
## slope (1/slope)                     -0.37986728  -0.03499134   0.30988461
## elev (1/elev)                       -0.04693585  -0.01277054   0.02139478
## area (square kilometers)           230.24434776 398.33132029 611.73252692
## τ[position] (days)                   4.56389649   7.82595173  13.41956827
## τ[velocity] (minutes)               40.67981657  43.18035547  45.83459945
## speed (kilometers/day)              11.89900607  13.45303449  15.00482017
## diffusion (square kilometers/day)    4.14577993   5.29369812   6.57977920

A suitability map

suitability_riemann <- ctmm::suitability(cilla_rsf_riemann, R=be, grid=crop(be[[1]], extent(cilla) * 2))
raster::plot(suitability_riemann)

Range distribution (includes the ranging behaviour)

my_grid <- crop(be[[1]], extent(cilla) * 2)
agde_cilla <- agde(CTMM = cilla_rsf_riemann, R = be, grid = my_grid)
plot(agde_cilla)

# Plot raster of range distribution
agde_raster <- ctmm::raster(agde_cilla, DF = "PMF")
plot(agde_raster)

Traditional RSF with downweighted Poisson regression

Functions to generate quadrature points and predict with the model

rsf_points <- function(x, UD, R = NULL, n = 1e5, k = 1e6, type = "Riemann", 
                       rmax =6*sqrt(UD@CTMM$sigma[1,1]),
                       interpolation = FALSE) {
  # Samples background points from a 2D normal distribution fitted to the relocation data, and extracts environmental
  # information from a raster object "R"
  # x: telemetry object
  # UD: UD object
  # R: raster* object
  # n: number of background points to sample
  # k: weight of presence points
  # rmax: maximum distance for Riemann-type integration
  # interpolation: do interpolation when sampling the grid
  # When type 0´= "MonteCarlo", importance sampling is done
  stopifnot(UD@CTMM$isotropic)
  stopifnot(type %in% c("Riemann", "MonteCarlo"))
  if (type == "Riemann") {
    quadrature_pts <- getValues(R)
    xy <- as.data.frame(xyFromCell(R, 1:ncell(R)))
    xy <- sf::st_as_sf(xy, coords = c("x", "y"), crs = crs(R))
    xy <- st_transform(xy, CRS(UD@info$projection))
    xy <- st_coordinates(xy)
    r <- sqrt(((xy[,1] - UD@CTMM$mu[1]))^2 + ((xy[,2] - UD@CTMM$mu[2]))^2)
    bg <- data.frame(case_ = 0,
                     x_ = xy[r<rmax,1], y_ = xy[r<rmax,2],  w_ = prod(res(R)), k_ = k)
    bg <- cbind(bg, quadrature_pts[r<rmax,])
    bg <- sf::st_as_sf(bg, coords = c("x_", "y_"), crs = CRS(UD@info$projection))
    xx <- data.frame(case_ = 1, x_ = x$x, y_ = x$y,
                     w_ = 1/k * UD$weights * mean(UD$DOF.area), k_ = k)
    xx <- sf::st_as_sf(xx, coords = c("x_", "y_"), crs = CRS(UD@info$projection))
    xx[names(R)] <- as.data.frame(raster::extract(R, st_transform(xx, crs(R)), method = ifelse(interpolation, "bilinear", "simple")))
    xx <- rbind(bg, xx)
  } else {
    quadrature_pts <- MASS::mvrnorm(n, mu = UD@CTMM$mu, Sigma = UD@CTMM$sigma)
    xx <- data.frame(case_ = 0, x_ = quadrature_pts[, 1], y_ = quadrature_pts[, 2], w_ = UD@CTMM$sigma[1,1]/n, k_ = k)
    xx <- rbind(xx, data.frame(case_ = 1, x_ = x$x, y_ = x$y,
                               w_ = 1/k * UD$weights * mean(UD$DOF.area), k_ = k 
    ))
    xx <- sf::st_as_sf(xx, coords = c("x_", "y_"), crs = CRS(UD@info$projection))
    xx[names(R)] <- as.data.frame(raster::extract(R, st_transform(xx, crs(R)), method = ifelse(interpolation, "bilinear", "simple")))
  }
  xy <- st_coordinates(xx)
  colnames(xy) <- c("x_", "y_")
  sd <- sqrt(UD@CTMM$sigma[1,1])
  xy[,1] <- (xy[,1] - UD@CTMM$mu[1])/sd
  xy[,2] <- (xy[,2] - UD@CTMM$mu[2])/sd
  xx <- cbind(xy, xx)
  xx
}

predict_rsf <- function(model, UD, object, include_avail = TRUE, data_crs = UD@info$projection) {
  if(include_avail) {
    xy <- as.data.frame(xyFromCell(object, 1:ncell(object)))
    xy <- sf::st_as_sf(xy, coords = c("x", "y"), crs = crs(object))
    xy <- st_transform(xy, data_crs)
    xy <- st_coordinates(xy)
    colnames(xy) <- c("x_", "y_")
    
    sd <- sqrt(UD@CTMM$sigma[1,1])
    xy[,1] <- (xy[,1] - UD@CTMM$mu[1])/sd
    xy[,2] <- (xy[,2] - UD@CTMM$mu[2])/sd
    
  } else {
    xy <- cbind("x_" = rep(0, ncell(object)), "y_" = rep(0, ncell(object)))
  }
  newdata <- as.data.frame(cbind(xy, getValues(object)))
  lambda <- as.numeric(predict(model, newdata, type = "link"))
  r <- raster(object, 1)
  r[] <- exp(lambda - max(lambda, na.rm = TRUE))
  r <- r / sum(getValues(r), na.rm = TRUE)
  r
}

A minimal “classic” example

Generate quadrature points (“background points”)

set.seed(2)
rsf_cilla_df <- rsf_points(cilla_telemetry, cilla_akde, buffalo_env, interpolation = TRUE)
rsf_cilla_df <- rsf_cilla_df[!is.na(rsf_cilla_df$slope),]

Fit a downweighted Poisson regression

m_rsf_cilla <- glm(case_*k_ ~ x_ + y_ + I(-(x_^2 + y_^2)/2) + elev + slope + var_NDVI, 
                  family = poisson(), data= rsf_cilla_df, weights = w_)

## Warning: glm.fit: fitted rates numerically 0 occurred

Summary of model and confidence intervals for parameter estimates

summary(m_rsf_cilla)

## 
## Call:
## glm(formula = case_ * k_ ~ x_ + y_ + I(-(x_^2 + y_^2)/2) + elev + 
##     slope + var_NDVI, family = poisson(), data = rsf_cilla_df, 
##     weights = w_)
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -10.32151    3.56117  -2.898  0.00375 ** 
## x_                    0.31218    0.36103   0.865  0.38721    
## y_                    0.48703    0.44872   1.085  0.27775    
## I(-(x_^2 + y_^2)/2)   1.14409    0.27988   4.088 4.35e-05 ***
## elev                 -0.01275    0.01742  -0.732  0.46426    
## slope                -0.03507    0.17592  -0.199  0.84201    
## var_NDVI            -72.71439   58.86105  -1.235  0.21670    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1125.8  on 289060  degrees of freedom
## Residual deviance: 1057.5  on 289054  degrees of freedom
## AIC: 1071.5
## 
## Number of Fisher Scoring iterations: 24

confint(m_rsf_cilla)

##                            2.5 %     97.5 %
## (Intercept)          -16.6341489 -2.9214401
## x_                    -0.3812531  1.0607780
## y_                    -0.3346315  1.4178533
## I(-(x_^2 + y_^2)/2)    0.6772977  1.7773232
## elev                  -0.0487675  0.0168988
## slope                 -0.4896424  0.2023228
## var_NDVI            -185.7287720 45.8497752

Map of suitability, including the home ranging behaviour

suitability_glm <- predict_rsf(m_rsf_cilla, cilla_akde, crop(buffalo_env, extent(cilla) * 2))
raster::plot(suitability_glm)

Map of suitability, without the home ranging behaviour

suitability_no_avail_glm <- predict_rsf(m_rsf_cilla, cilla_akde, crop(buffalo_env, extent(cilla) * 2), include_avail = FALSE)
raster::plot(suitability_no_avail_glm)

Comparison of estimates from the two approaches (ctmm and “classic” via glm)

vars <- c("elev", "slope", "var_NDVI")

# Estimates are not identical but similar
cilla_rsf_riemann$beta[vars]

##         elev        slope     var_NDVI 
##  -0.01277054  -0.03499134 -72.68967041

coef(m_rsf_cilla)[vars]

##         elev        slope     var_NDVI 
##  -0.01274543  -0.03506537 -72.71439217

# Confidence intervals are of similar width, but location is different

ci_glm <- confint(m_rsf_cilla)
ci_glm[vars, ]

##                 2.5 %     97.5 %
## elev       -0.0487675  0.0168988
## slope      -0.4896424  0.2023228
## var_NDVI -185.7287720 45.8497752

summary(cilla_rsf_riemann)$CI[sapply(vars, function(x) grep(x, rownames(summary(cilla_rsf_riemann)$CI))), ]

##                                 low          est        high
## elev (1/elev)           -0.04693585  -0.01277054  0.02139478
## slope (1/slope)         -0.37986728  -0.03499134  0.30988461
## var_NDVI (1/var_NDVI) -188.04790701 -72.68967041 42.66856618

Create maps of suitability based on parameter estimates

M <- getValues(buffalo_env)
eta_glm <- M[,vars] %*% coef(m_rsf_cilla)[vars]
eta_rsf <- M[,vars] %*% cilla_rsf_riemann$beta[vars]

suit_r_glm <- raster(buffalo_env, 1)
suit_r_glm[] <- exp(eta_glm)
suit_r_rsf <- raster(buffalo_env, 1)
suit_r_rsf[] <- exp(eta_rsf)

suit_raster <- raster::stack(suit_r_rsf, suit_r_glm)
names(suit_raster) <- c("ctmm", "classic")
raster::plot(suit_raster)

raster::plot(suit_raster, ext = extent(cilla) * 2)

Multiple animals

• with rsf.fit: you can use the mean function on a list of rsf.fit objects
• “classic” approach: use glmmTMB and a mixed-effects Poisson regression: Muff, Signer & Fieberg (2020) J Anim Ecol 89: 80-92.